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Abstract 

The interplay between higher orders of the perturbative QCD (pQCD) expansion and higher- 
twist contributions in the analysis of recent Jefferson Lab data on the lowest moment of the 
spin-dependent proton Ti(Q 2 ) at 0.05 < Q 2 < 3GeV 2 is studied. We demonstrate that the 
values of the higher-twist coefficients f£'r extracted from the data by using the singularity-free 
analytic perturbation theory provide a better convergence of the higher-twist series than with the 
standard perturbative QCD. From the high-precision proton data, we extract the value of the 
singlet axial charge ao(l GeV 2 ) = 0.33 ± 0.05. We observe a slow Q 2 dependence of fitted values of 
the twist coefficient /X4 and ao when going to lower energy scales, which can be explained by the 
renormalization group evolution of /^(Q 2 ) and ao(Q 2 ). As the main result, a good quantitative 
description of all the Jefferson Lab data sets down to Q ~ 350 MeV is achieved. 
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I. INTRODUCTION 



The spin structure of the nucleon remains the essential problem of nonperturbative QCD 
and hadronic physics. One of its most significant manifestations is the so-called spin crisis 
or spin puzzle related to the surprisingly small fraction of proton polarization carried by 
quarks [l|, 0]. This problem attracted attention to the peculiarities of the underlying QCD 
description of the nucleon spin, in particular, to the role of the gluonic anomaly (see [l|, Q 
and references therein). The natural physical interpretation of these effects was the gluon 
(circular) polarization, while the experimental indications of its smallness may also point to 
a possible manifestation of the anomaly via the strangeness polarization j^]. The key point 
is its consideration as a kind of heavy-quarks polarization [BJ due to the multiscale 0] picture 
of the nucleon exploring the fact that strange quark mass is much (as the squared ratios 
matter) smaller than the nucleon one and, in turn, larger than higher-twist parameters. 

Higher-twist parameters (known also as the color polarizabilities) are important ingredi- 
ents of the nucleon spin structure. Their extraction from experimental studies is relatively 
complicated as they are most pronounced at low momentum transfer Q. Although in this 
region very accurate Jefferson Lab (JLab) data are now available, higher-twist contributions 
are shadowed by Landau singularities of QCD coupling. As was shown in Ref. jsj, this prob- 
lem may be solved by the use of singularity-free couplings which allowed a quite accurate 
extraction of higher twist (HT) and a fairly good description of data down to rather low Q. 
The object of investigation in j6[ was the difference of the lowest moments rf' n of proton 
and neutron structure functions g±, which corresponds to the renowned Bjorken sum rule 
(BSR) [tJ. At finite Q 2 the moments T^ n are modified by higher order radiative corrections 
and higher- twist power corrections, as dictated by the operator product expansion (OPE). 
Such generalized (Q 2 -dependent) BSR became a convenient and renowned target ground for 
testing different possibilities of combining both the perturbative and nonperturbative QCD 
contributions in the low-energy domain (see, for example, Refs. 0,0]). 

The global higher twist analysis of the data on the spin-dependent proton structure 
function g\ at relatively large 1 < Q 2 < 30 GeV 2 , was performed in Ref. [hJ. While the 1/Q 2 
term in the OPE works at relatively high scales Q 2 > 1 GeV 2 , higher- twist power corrections 
1/Q 4 , l/Q 6 , etc., start to play a significant role at lower scales, where the influence of the 
ghost singularities in the coefficient functions within the standard perturbation theory (PT) 
becomes more noticeable. It affects the results of extraction of the higher twists from 
the precise experimental data leading to unstable OPE series and huge error bars [6]. It 
seems natural that the weakening or elimination of the unphysical singularities of the QCD 
coupling would allow shifting the perturbative QCD (pQCD) frontier to a lower energy scale 
and getting more exact information about the nonperturbative part of the process described 
by the higher-twist series. 

As was shown in Ref. the situation becomes better if one uses for a running coupling 
a more precise iterative solution of the renormalisation group (RG) equation in the form of 
the so-called denominator representation [l2| instead of the Particle Data Group loop 1/L 
expansion [l3|, especially at the two- loop level. In this investigation, to avoid completely 
the unphysical singularities at Q = Aqcd ~ 400 MeV we deal with the ghost-free analytic 
perturbation theory (APT) [jjj (for a review on APT concepts and algorithms, see also 
Ref. [3]), which recently proved to be an intriguing candidate for a quantitative description 
of light quarkonia spectra within the Bethe-Salpeter approa ch [16[| . and the so-called glueball- 
freezing model proposed recently by Yu. A. Simonov in Ref. [17| (below, SGF model) to avoid 
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the renormalon |18l| ambiguity in QCD. Other versions of frozen a s models were developed 
earlier in Ref. [191 ] . As it will seen below that APT and SGF approaches predict very 
close couplings at Q > Aqcd, whereas they have different infrared-stable points at Q = 0. 
Consequently, as it was shown in Ref. [6[ , these models lead to very close perturbative parts 
of the Bjorken sum rf~™ rt . The higher-twist contributions turned out to be very close, too. 
Here, we would like to discuss this point in more detail. 

In the current paper we study the interplay between higher orders of the pQCD expansion 
and higher-twist contributions using the recent JLab data on the lowest moments of the 
spin-dependent proton and neutron structure functions Y^ n {Q 2 ) and Y\~ n {Q 2 ) in the range 



0.05 < Q < 3 GeV [11]. Thus, we extend and generalize the analysis started in Ref. [6| by 
considering also the singlet channel involving the Y p { n (Q 2 ) for the proton (providing the most 
accurate data) and the neutron structure functions separately. This allows, in particular, 
determining the singlet axial charge a coming into both Y^ n (Q 2 ) moments, which in the 
quark-parton model is identified with the total spin carried by quarks in the proton. For this 
purpose, we perform the global analysis of the JLab precise low-energy data on Y^(Q 2 ) [2(| 
using the advantages of the APT and SGF model, and extract the singlet axial charge a , as 
well as the coefficient fj,^' n of the 1/Q 2 subleading twist-4 term, which contains information 
on quark-gluon correlations in nucleons. 

The paper is organized as follows. In Sec. 2, the lowest moments analysis for the polarized 
structure functions g v { n in the framework of the conventional PT approach is performed. In 
Sec. 3, we dwell briefly on the APT, its ideas and the results of its application to Y^ n {Q 2 ). 
In Sec. 4, we apply the formalism to the analysis of the low-energy data on the first 
moments Y p,n (Q 2 ) and compare the results with the results of other researchers concerning 
the singlet axial constant and gluon polarization Ag at low Q 2 < 1 GeV 2 . Section 5 
contains discussion and some concluding remarks. 



II. SPIN SUM RULES IN CONVENTIONAL PT 
A. First moments of spin structure functions g^' n 

The lowest moments of spin- dependent proton and neutron structure functions g v { n are 
defined as follows: 

i 



rr(Q 2 )= / dxgr(x,Q 2 ), (2.1) 
Jo 

with x = Q 2 j2Mv, the energy transfer z/, and the nucleon mass M. The upper limit includes 
the proton/neutron elastic contribution at x = 1. This contribution becomes essential if the 
OPE is used to study the evolution of the integral in the moderate and low momentum 



transfer region Q < 1 GeV [2JJ . It is of special interest to analyze data with the elastic 
contribution excluded, since the low-Q 2 behavior of "inelastic" contributions to their nons- 
inglet combination Y^~ n (Q 2 ), i.e. BSR, is constrained by the Gerasimov-Drell-Hearn (GDH) 
sum rule j22| , and one may investigate its continuation to a low scale [s| . So below we study 
inelastic contributions Y v ^ el ^Q 2 ) using the corresponding low-energy JLab data jloj. Note 
that the influence of the "elastic" contribution is noticeable starting from the higher-twist 
~ fiQ term which is natural due to a decrease of the elastic contribution with growing Q 2 [f| . 

At large Q 2 the moments Y^ ,n (Q 2 ) are given by the OPE series in powers of 1/Q 2 with 
the expansion coefficients related to nucleon matrix elements of operators of a definite twist 
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(defined as the dimension minus the spin of the operator), and coefficient functions in the 
form of pQCD series in a™ (see, e.g., Ref. 23)). In the limit Q 2 3> M 2 the moments are 
dominated by the leading twist contribution, $ n {Q 2 )-, which is given in terms of matrix 
elements of the twist-2 axial vector current, ip r )^ r )^ip. This can be decomposed into fla- 
vor singlet and nonsinglet contributions. The total expression for the perturbative part of 
r^ n (Q 2 ) including the HT contributions reads 



rr(Q 2 ) = ^ 
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where Es and E^s are the singlet and nonsinglet Wilson coefficients, respectively, calculated 
as series in powers of a s |24|. These coefficient functions for n/ = 3 active flavors in the MS 
scheme are 
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The triplet and octet axial charges as = gA— 1.267 ±0.004 [13| and a§ = 0.585 ±0.025 [25], 
respectively, are extracted from weak decay matrix elements and are known from /3-decay 
measurements. As for the singlet axial charge ao, it is convenient to work with its RG 
invariant definition in the MS scheme a 1 ™ = a Q (Q 2 = oo), in which all the Q 2 dependence 
is factorized into the definition of the Wilson coefficient Es(Q 2 ). 

In contrast to the proton and neutron spin sum rules (SSRs), the singlet and octet 
contributions are canceled out, giving rise to more fundamental BSR 
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ri- n m = 9 fE NS (Q 2 )+ Y^^M 
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which is analyzed here along with the proton SSR in more detail than in Ref. [6(. The first 
nonleading twist term [26| can be expressed (27| 
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in terms of the color polarizability fi- 

The RG Q 2 evolution of the axial singlet charge ao(Q 2 ) and nonsinglet higher- twist 



Ft n (Q 2 ) is 
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In the NLO we may write 

a (Q 2 ) ~ a (Q 2 ) [l + A^Q 2 ) + 0(a 



(2.6) 
(2.7) 

(2.8) 
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As a first step of our analysis, in Eq. ( 12. 2 j) we will neglect the weak dependence of //^ n on 
logC} 2 - Note that the evolution of the higher-twist terms //e,8,... i n Eq. (12. 2p is still unknown. 
As a next step we discuss the possible influence of the H^{Q 2 ) evolution on our results. 
The Q 2 evolution of the proton higher- twist term n\{Q 2 ) is assumed to be the same as the 
evolution of the nonsinglet twist tfl~ n (Q 2 ). This may be justified by the relative smallness 
of the singlet higher-twist term. 



TABLE I: Current NLO fit results for the axial singlet charge uq. 



Reference 


LSS [30] 


DSSV [31] 


AAC [32] 


HERMES [29] 


COMPASS [28] 


Ql GeV 2 


1.0 


10.0 


4.0 


5.0 


3.0 


a 


0.24 ± 0.07 


0.24 


0.25 ±0.05 


0.32 ± 0.04 


0.35 ±0.06 



Let us discuss current results for the nucleon spin structure and higher twists. In Table [H 
we list the fit results for the axial singlet charge ao from the literature including all global 
NLO PT a naly ses and the recent results obtained directly from deuteron data on Tf by 
COMPASS [28] and HERMES @- The global fit results for a are somewhat lower than 



that from the deuteron data. It was mentioned in the most recent review [2] that the reason 
for such a discrepancy is not completely understood. Further, we analyze this issue in more 
detail. 



TABLE II: Current NLO fit results for the highest-twist term /U4/M 2 . The uncertainties are statistical 
only. 



Target 


Proton [36] 


Neutron [37] 


p-n [38] 


p-n [35] 


p - n [6] 


Q 2 , GeV 2 


0.6 - 10.0 


0.5 - 10.0 


0.5 - 10.0 


0.66 - 10.0 


0.12 - 3.0 


//4/M 2 


-0.065 ±0.012 


0.019 ±0.002 


-0.06 ± 0.02 


-0.04 ±0.01 


-0.048 ± 0.002 



A detailed higher-twist analysis based on the combined SLAC and JLab data [on proton, 



neutron Ti' n (Q 2 ) [33| and nonsinglet T^^C} 2 ) moments 35]] was performed in Refs. 351-138 
In Table [TT1 we show the current results for the twist-4 coefficient /Z4/M 2 at Q 2 = 1 GeV 2 
extracted from Ff' n data. As we have seen from our previous analysis [6|, a satisfactory 
description of the low-energy JLab data on the Bjorken sum rule down to Q m in 
350 MeV can be achieved by using APT and taking into account only three higher-twist 
terms fj^ ~ 6 n 8 . Including only the twist-4 term fJ^~ n /M 2 , this method allowed us to get its value 
with noticeably higher accuracy than in the standard PT approach, shifting the applicability 
of the pQCD expansion down to Q 2 min = 0.47 GeV 2 . The higher-twist analysis of the most 
recent precise JLab experimental data on the proton spin sum rule [20|] has not been carried 
out yet in the literature. This gives us a reasonable motivation for a detailed data analysis 
and studying the higher-twist effects at low-energy scale both in the standard PT, APT and 
"infrared-frozen" a s approaches. 
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B. The running coupling 



The infrared behavior of the strong coupling is crucial for the extraction of the nonper- 
turbative information from the low-energy data. Within the pQCD, the a s coupling can be 
found by a solution of the RG equation 



da s 



-/3 a 2 (l + bia s + b 2 a 2 s + ... 



where L = ln(Q 2 /A 2 ) and b^ = Pk/Po- The standard PT running coupling a s is usually 
taken in the form [see, for example, Eq. (6) in the recent review 39] or Eq. (9.5) in the PDG 
review 



13| ] expanded in a series over In L/L, i.e 
1 



fexlnL 1 



[6 2 (ln 2 L-lnL-l) + 6 2 ] . 



(2.9) 



Here, the 1/L 2 term corresponds to the 2-loop contribution and the 1/L 3 term is usually 
referred to as "the 3-loop one." Actually, the pieces of genuine 2-loop contribution pro- 
portional to bi are entangled with the higher-loop ones. This defect is absent in the more 
compact denominator representation [l2{|, which at 2, 3-loop levels has the following forms: 
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which, being generic for the PDG expression (I2.9p . are closer to the corresponding iterative 
RG solutions and, hence, more precise. Advantages of formulas (I2.10p in__the higher-twist 
analysis of the Bjorken sum rule were demonstrated in our previous work 




FIG. 1: The NLO running coupling a s in different approaches. 



In Fig. [U we compare the behavior of the two-loop running coupling a s at low Q 2 scales 
in different approaches. The long-dashed line is the exact two-loop PT result, the dotted 
line is the denominator representation (I2.10p (referred to as "Denom" below), and the short- 
dashed line is the PDG expression (12. 9p . As one can see from this figure, the NLO Denom 
coupling is much closer to the corresponding numerical RG solution than the 1 / L-expanded 
PDG expression. 
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In Fig. [TJ we also show two models of the infrared-stable running coupling. One of them 
is the Simonov "glueball-freezing model" (SGF-model) [l7j], represented by the dash-dotted 
line, with the 1/L-type loop expansion for the "infrared-frozen" coupling similar to PDG 

a B (Q 2 ) = ^(L), L = ln(?-±M), (2.11) 

where the two-loop as is taken in the form of the first two terms in Eq. ( 12.91) with logarithm 
modified by a "glueball mass" M ~ 1 GeV. Note, the usual PT expansion in powers of 
as in the coefficient functions ( 12. 3 p and (12. 4p is adopted. The solid line corresponds to the 
second model of the infrared-stable coupling - the APT running coupling, which will be 
discussed in detail below in the next section. 

As one can see from Fig. [TJ the SGF and APT couplings are very similar in the low-energy 
domain Aqcd < Q 1 GeV though their infrared limits are different. Also, a comparison 
of APT and PT couplings over a wide range of Q 2 , 1 < Q 2 < 10 4 GeV 2 , can be found in 
Ref. 0. 

Note, we extract values of Aqcd corresponding to different models of the running cou- 
pling, by evolution from the world experimental data on a s (M|) as a normalization point 
in each particular order of PT. 



C. Stability and duality 

In the following, when calculating the observables in any particular order of perturbation 
theory, we will employ the prescription for the coefficient functions in the infrared region, 
where the order of the power a s series in the coefficient functions is matched with the loop 
order in a s itself. For example, for the nonsinglet coefficient function in the Bjorken sum 



rule, we write consequently (for details, see Ref. |29| ) 



E L N °s = 1, E™ = 1 - ^— , < L ° = 1 - ^ 3.558 (^—) , . . . (2.12) 

We see that the leading singular behavior in the coefficient function ~ In™ L/L m when L — >■ 
comes from the highest power of a s . So in the infrared domain the influence of singularities 
gets stronger in higher orders of perturbation theory that may affect the data analysis below 
1 GeV 2 . This fact explains our observation made in Ref. 0], where we showed that the higher 
PT orders yield a worse description of the BSR data in comparison with the leading order. 
We observe a similar picture for the precise JLab data on T^(Q 2 ) [20] probably implying the 
asymptotic character of the series in powers of a s (see Fig. [2]). 

The corresponding fit results for HT terms, extracted in different orders of PT, are listed 
in Table IIII1 We see that with raising the loop order the values of fJ^ g terms increase, 
whereas //g decreases, yielding a "swap" between the higher orders of PT and HT terms. 
Such a "swap" between PT and HT terms (decreasing HT term by including more terms of 



PT and using resummation of PT series) was previously observed in Refs. [41|, |42J. A similar 
situation holds when fitting T^(Q 2 ) data over the fixed range 0.8 GeV < Q < 2.0 GeV , where 
it is sufficient to take into account only one twist term /i 4 . 

In Fig. [31 we show fits of BSR data (left panel) and proton SSR data (right panel) in 
different orders of perturbation theory taking only into account the fi4 term. One can see 
there that the higher-loop contributions are effectively "absorbed" into the value of /14 which 
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FIG. 2: Best fits of JLab and SLAC data on BSR r?""(Q 2 ) (left panel) and proton SSR T\(Q 2 ) (right 
panel) calculated at various loop orders. 

TABLE III: Dependence of the best fit results of BSR T\~ n (Q 2 ) and proton SSR T\{Q 2 ) data (elastic 
contribution excluded) on the order of perturbation theory [NLO and NNLO Denom couplings (12.101) 
are used]. The corresponding fit curves are shown in Fig. [2J The minimal borders of fitting domains 
in Q 2 are settled from the ad hoc restriction x 2 ^ 1 ar, d monotonous behavior of the resulting fitted 
curves. 



Target 


Method 


QLn, GeV 2 


inv 

"o 


^/M 2 


W/M 4 


/i 8 /M 6 




LO 


0.121 


0.29(2) 


-0.089(3) 


0.016(1) 


-0.0010(1) 


proton 


NLO 


0.17 


0.38(2) 


-0.070(5) 


0.010(2) 


0.0004(3) 




NNLO 


0.38 


0.37(5) 


-0.034(19) 


-0.025(20) 


0.017(6) 




LO 


0.17 




-0.126(5) 


0.037(3) 


-0.004(1) 


p - n 


NLO 


0.17 




-0.076(5) 


0.019(3) 


-0.001(1) 




NNLO 


0.38 




-0.026(11) 


-0.035(15) 


0.026(5) 



decreases in magnitude with increasing loop order while all the fitting curves are very close 
to each other. This observation reveals a kind of "duality" between the perturbative a s 
series and nonperturbative 1 / Q 2 series. A similar phenomenon was observed before for the 
structure function F-* in Refs. 



43. 44 



This also means the appearance of a new aspect of quark hadron duality, the latter being 
the necessary ingredient of all the QCD applications in the low-energy domain. Usually, it 
is assumed [45| that the perturbative effects are less important there than the power ones 
due to a nontrivial structure in the QCD vacuum. 

In our case, the PT corrections essentially enter into the game, so that the pQCD higher 
order terms are relevant in the domain where the concepts of traditional hadronic physics 
are usually applied. 

The interplay between partonic and hadronic degrees of freedom in the description of 
GDH SR and BSR may also be observed in the surprising similarity between the results of 
"resonance" {46[ and "parton" j9[ approaches. 



S 



rr(Q) 



r/(Q) 




JLab Hall B (CLAS EG 1b) [2006] 
(stat. errors only) 

JLab Halls A,B E94010/EG1a [2002] 
JLab Hall B (CLAS EG1a) [2003] 
S LAC E 143 [2000] 



LO Bjorken: = 


-0.084(4} 


NLO "Denom": n"" 


= -0.045(3) 


NNLO "Denom": ^ 


"" = -0.019(4) 




LO proton SSR: ^ 


= -0.069(2), a o 


= 0.24(2) 


NLO "Denom": (i a P 


= -0.057(2), a = 


0.35(2} 


NNLO "Denom": (i 


p = -0.038(2), a D 


= 0-35(2) 



Q, GeV 



Q, GeV 



FIG. 3: Best fits of JLab and SLAC data on BSR Tf 
panel) calculated in various loop orders with fixed Qmin 



\Q 2 ) (left panel) and proton SSR T\{Q 2 ) (right 
= 0.8 GeV. 



TABLE IV: Dependence of the best (3 + l)-parametric fit results of I^(Q 2 ) data (elastic contribution 
excluded) on A„ /= 3 in NLO Denom PT. 



Aqcd, MeV 


QL. GeV 2 


„inv 
«0 




/VM 4 


fi 8 /M 6 


300 


0.14 


0.40(2) 


-0.077(3) 


0.014(1) 


-0.0005(2) 


400 


0.24 


0.39(3) 


-0.064(8) 


0.006(5) 


0.002(1) 


500 


0.35 


0.34(4) 


-0.028(13) 


-0.033(11) 


0.019(3) 



One may ask to what extent these results are affected by the unphysical singularities 
when approaching Q ~ Aqcd in the PT series for r^'p T . Their influence becomes essential 
at Q < 1 GeV where the HT terms start to play an important role. The minimal border 
of the fitting domain Q m i n is tightly connected with the value of Aq C d', i-e. it is a scale, 
below which the influence of the ghost singularities becomes too strong and destroys the 
fit. To see how the Q 2 min scale and fit results for the \i terms change with varying Aqcd, 
we have performed three different NLO fits with Aqcd = 300, 400, 500 MeV (see Table HVJ). 
It turns out that the term /z 4 is quite sensitive to the Landau singularity position, and its 
value noticeably increases with increasing Aqcd- The APT and "soft-frozen" models are 
free of such a problem, thus providing a reliable tool of investigating the behavior of HT 
terms extracted directly from the low-energy data j6[. This provides a motivation for the 
analysis performed in the next section. 



III. MOMENTS If n (Q 2 ) IN ANALYTIC PERTURBATION THEORY 

The moments of the structure functions are analytic functions in the complex Q 2 plane 
with a cut along the negative real axis, as was demonstrated in Ref. 1471 (see also Ref. |48fl). 
On the other hand, the standard PT approach does not support these analytic properties. 
The influence of requiring these properties to hold in the DIS description was studied previ- 
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ously by Igor Solovtsov and coauthors in Refs. [40|, |49j. Here we continue this investigation 
by applying the APT method, which gives the possibility of combining the RG resummation 
with correct analytic properties of the QCD corrections, to the low-energy data on nucleon 
spin sum rules Y^ n {Q 2 ). 

In the framework of the analytic approach we can write the expression for T^ n (Q 2 ) in 
the form 



r p,n 
1 I, APT 



(Q 2 



1 

12 



±a 3 + la s yr s T (Q 2 ) + \< v Ef T m 



oo APT;p,n 



lh 



(Q 2 



8=2 



Q 



2i-2 



(3.1) 



which is analogous to one in the standard PT (|2.2p . The corresponding NNLO APT modi- 
fication of the singlet and nonsinglet coefficient functions is 



E 



J NS 
APT 



(Q 2 



1 -0.318 Af\Q 2 ) -0.361 A ( i\Q 2 ) - ... 
1 -0.318 Af{Q 2 ) -0.111 A ( i\Q 2 ) - ... 



(3.2) 
(3.3) 



(3) 

where A k is the analyticized fcth power of 3-loop PT coupling in the Euclidean domain 



At\Q 2 
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o~ + Q 2 



n = 3 . 



(3.4) 



In the one-loop case, the APT Euclidean functions are simple enough [14 



A? {I) 
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Jo 
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L + A^ 
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Wo 
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L 2 



Q 2 
Q 2 A 2 



In 



A 2 



(3.5) 



(Q 2 - A 2 y 



A (l) 



1 dA 



(i) 



k Po dL 



i.e. the higher functions Ak are related to the lower ones recursively by differentiating. 
Analogous two- and three-loop level expressions involve the special Lambert function and 
are more intricate, and they can be found in Refs. [5(J 51]. It should be stressed that the APT 
couplings are stable with respect to different loop orders at low-energy scales Q 2 < 1 GeV 2 
[l5| . This feature is absent in the standard PT approach, as reflected in Fig. [2j 

Meanwhile, even for the three-loop APT case, there exists a possibility to employ the 
effective log approach proposed by Igor Solovtsov and one of the authors in Ref . [52j • In the 
present context, in the region Q < 5 GeV one may use simple model one-loop expressions 
(13. 5p with some effective logarithm L* : 



A?l 3 (L) -+ AT£ = A\%(L*) , — 2 HQ/ A 



i(i) 



«), A^0.50A( 3 ). 



(3.6) 



Thus, instead of the exact three-loop expressions for the APT functions, in Eq. f)3.3p one 
can use the one-loop expressions 03.51) with the effective A parameter A mod = A^ whose 
value is given by the last relation (13. 6p . This model was successfully applied for higher-twist 
analysis of low-energy data on BSR in our previous work 0], and also in the T decay analysis 
in Ref. @. 

The maximal errors of the model (13.61) for the first and the second functions are 

A n =3 , which seem to be sufficiently 



5A? od /Ar ~ 4% and 5 AT /A 
accurate. Indeed, as far as A\{Q 



mod 



8% at Q 
400 MeV) 



0.532 and ^ 2 (400MeV) = 0.118, 
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FIG. 4: Evolution of do (Q 2 ) normalized at Qo = FIG. 5: Evolution of fj% ™(Q 2 ) normalized at 
lGeV 2 . Q 2 = lGeV 2 . 



the total error in T^ APT is mainly determined by the first term, being of the order 
ST^/T^ ~ 5A™ od /n ~ 1 % , i.e., less than the data uncertainty. 

In order to take into account the one-loop Q 2 evolution of the axial singlet charge a (Q 2 ), 
we use expression (12.91) substituting the one- loop analytic coupling A[ (L). The contribution 
of the ~ Ai term to ao(Q 2 ) at, for example, Q 2 = 0.1 GeV 2 with normalization point at 
Ql = 1 GeV 2 is A^O.lGeV 2 ) ~ 0.11; i.e. the evolution contributes about 10% when one 
shifts the pQCD border down to Aqcd (see Fig. HJ). 

For the evolution of the twist-4 term ^^{Q 2 ) (12. 7p . we have to "analyticize" the fractional 
power (a s ) u . For this purpose we apply the fractional APT approach developed in Ref. 54 
At the one-loop level in the Euclidean domain we have 



1 

J? 



F(e 



(3.7) 



Here F(z, v) is the Lerch transcendent function. In this case, the evolution of the nonsinglet 
twist-4 term in BSR reads 



p—n I r\i 
A*4, APT W 



p—n 
^4,APT 



(qi: 



A {1) 



(Q 2 



A v (Qq 



V 



32 
81' 



(3.8) 



The corresponding evolution is shown in Fig. |5j As follows from this figure, the evolution 
from 1 GeV to Aqcd increases the absolute value of 1^2pt by about 20 %. 



IV. NUMERICAL RESULTS 



A. Nonsinglet case: the Bjorken sum rule 



In Fig. El we show best fits of the combined data set for the BSR function T I l~ n (Q 2 ) 
with NLO Denom (solid lines) and PDG (dashed lines) couplings and NNLO APT (dash- 
dotted lines) at fixed Aq C d value corresponding to the world average. We also show here the 



11 



0,20 



0,18 

0,16 -|:| 
0,14 
0,12 
0,10 
0,08 
0,06 
0,04 
0,02 


-0,02 



r 



P-n 



NLO "Denom" 

NLO PDG 

- - N 2 LO APT 

SGF-model without twists 

N 2 LO APT without twists 

GDH slope 

GDH + matching with 

APT fit (a la Ref. [9]) 



,4 



i 




• JLab Hall B (CLAS EG1 b) [2006] 
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♦ JLab Hall B (CLAS EG1a) [2003] 
■ SLACE143[2000] 
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FIG. 6: Best 1,2,3-parametric fits of the JLab and SLAC data on Bjorken SR calculated with different 
models of running coupling. 



TABLE V: Combined fit results of BSR for the HT terms in APT, the SGF model and the standard 
PT approach. 



Method 


O 2 GeV 2 


/x 4 /M 2 


W/M 4 


/i 8 /M 6 




0.50 


-0.043(3) 








NLO PDG 


0.30 


-0.074(3) 


0.026(7) 







0.27 


-0.049(4) 


-0.010(3) 


0.010(1) 




0.47 


-0.049(3) 








NLO Denom 


0.17 


-0.069(4) 


0.014(1) 







0.17 


-0.065(7) 


0.011(3) 


0.0003(7) 




0.47 


-0.061(3) 








NLO SGF 


0.19 


-0.073(3) 


0.010(3) 







0.10 


-0.077(4) 


0.014(5) 


-0.0008(3) 




0.47 


-0.055(3) 








NNLO APT 


0.17 


-0.062(4) 


0.008(2) 





no evolution 


0.10 


-0.068(4) 


0.010(3) 


-0.0007(3) 




0.47 


-0.051(3) 








NNLO APT 


0.17 


-0.056(4) 


0.0087(4) 





with evolution 


0.10 


-0.058(4) 


0.0114(6) 


-0.0005(8) 



pQCD part of the BSR at different values of Ancn = 300, 400, 500 MeV calculated within 
APT (short-dashed lines) and the SGF model [17| at different values of the glueball mass 
M = 1.2, 1.0, 0.8 GeV (with A = 360 MeV) (dotted lines). 

The corresponding numerical results are given in Table E] As we have seen before in 
Fig. [TJ the behavior of SGF and APT couplings is very similar in the low-energy domain 
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Aqcd < Q ;$ 1 GeV. As a result, the corresponding perturbative parts of BSR in Fig. [6] 
and results for higher-twist terms in Table IVl turn out to be close, too. Our fits in APT and 
the SGF model give the HT values indicating a better convergence of the OPE series due to 
decreasing magnitudes and alternating signs of consecutive terms, in contrast to the usual 
PT fit results. 

As is seen from Table [Vj there is some sensitivity of fitted values of /i 4 with respect to 
Qmin variations; namely, it increases in magnitude when one incorporates into the fit the data 
points at lower energies. This property of the fit may be treated as the slow (logarithmic) 
evolution n^(Q 2 ) with Q 2 which becomes more noticeable at broader fitting ranges in Q 2 , 
as discussed above. So for completeness we included in Table IVl APT fits for /i 4 (<5o) taking 
into account their RG evolution with Qq — 1 GeV as a normalization point. We see that the 
fit results become more stable with respect to Q m %n variations. 

However, there is still a problem with how to treat the evolution of higher-twist terms 
^6,8,.. {Q 2 ) which again may turn out to be important when one goes to lower Q 2 , since the 
fit becomes more sensitive to very small variations of /i6,8,.. with Q 2 . 

Note that the APT functions Ak contain the (Q 2 )~ h power contributions which effectively 
change the fitted values of \i terms. In particular, subtracting an extra (Q 2 ) -1 term induced 
by the APT series 

rCxVW 2 ) - f + / ( + M- + o (jL) 

6 vih(qvaw )/ Q \Q J 

0.058 (with evolution) from Table IV] we 



0.042, A^ ~ 0.18 GeV (4.1) 

that nicely correlates with the result in Ref. [35]: fi^~ n /M 2 ~ —0.045. This demonstrates 
the concert of the APT analysis with the usual PT one for the BSR data at Q 2 > 1 GeV 2 . 

We do not take into account RG evolution in /z 4 for the standard PT calculations since 
the only effect of that would be the enhancement of the Landau singularities by extra 
divergencies at Aq C d (see Fig. E]), whereas at higher Q 2 ~ 1 GeV 2 the evolution is negligible 
with respect to other uncertainties. In ghost-free models, however, the evolution gives a 
noticeable effect at low Q ~ Aqcd- Note that our previous result in Ref. j6[, obtained 
without taking into account the RG evolution, turned out to be slightly larger than (14. ip 
H4~ n /M 2 ~ —0.048 which is very close to the corresponding value obtained with the most 
precise Denom PT coupling and is shown in Table |V] 

B. Singlet case: spin sum rules T^' n and nucleon spin structure 

Turn now to the three-loop APT part of the proton moment apt(Q 2 )- value is quite 
stable with respect to small variations of A, in contrast to the huge instability of PT : it 
changes now by about 2% - 3% within the interval A^ 3 ) = 300 - 500 MeV . The same was 
previously observed for the Bjorken function T P ^ PT {Q 2 ) in Ref. [6j. Because of this fact the 
low-Q 2 data on T^(Q 2 ) cannot be used for determination of A in the APT approach. 



with x = 0.43 and using the value ^^ PT /M = — 
finally get 

^apt + *Aff/ ^-(1 GeV 2 ) 
M 2 M 1 
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Extending the analysis of Ref. |49J to lower Q 2 scales, we estimated the relative size of 
APT contributions to T^Q 2 ). It turned out that the third term ~ A3 contributes no more 
than 5% to the sum, thus supporting the practical convergence of the APT series. 

TABLE VI: Sensitivity of the best APT fit results of proton T\{Q 2 ) data (elastic contribution excluded) 
to A n =3 variations. The minimal fitting border is Q 2 min = 0.12 GeV 2 . 



Aqcd, MeV 


"0 


/i 4 /M 2 


W>/M 4 


/i 8 /Af 6 


300 


0.43(3) 


-0.082(4) 


0.015(9) 


-0.0009(5) 


400 


0.45(3) 


-0.081(4) 


0.015(9) 


-0.0009(5) 


500 


0.47(3) 


-0.080(4) 


0.014(9) 


-0.0009(5) 



To see how the numerical fit results are sensitive to A 



(n f =3) 



in APT, we fulfilled four 



different fits of the proton Y\(Q 2 ) data with Aqcd = 300, 400, 500 MeV as we did before in 
the standard PT. The results of these fits are shown in Table I VI I Comparing these results 
with the data from Table IIV} we see that the corresponding results in the standard PT are 
much more sensitive to A variations than ones in APT. 



0.16 
0.14- 

0.12- 
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0.04- 
0.02- 
■ 
-0.02 • 
-0.04 ■ 



r/(Q) 



NLO PDG 
NLO "Denom" 
NNLO APT 

NNLO APT without twists 
NLO GF-model without twists 
GDH Q 2 slope 
GDH + matching 
with APT fit (a la Ref. [9]) 



■W6,8 ^4,6,8 // ^4,6 




• JLabCLAS EG1b 
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■ SLAC E143 
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0.3 



0.9 1.0 
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FIG. 7: Best (l,2,3+l)-parametric fits of the JLab and SLAC data on L^ (elastic contribution excluded). 



In Fig. we show best fits of the combined data set for the function T^Q 2 ) (the data 
uncertainties are statistical only) in the standard PT (PDG and Denom versions) and the 
APT approaches. We have also shown the perturbative parts of Y\{Q 2 ) calculated in APT 
and the SGF model. They are close to each other down to Q ~ A, similar to the BSR 
analysis in the previous subsection. A similar observation was made in the analysis of the 
small x spin averaged structure functions in Ref. [34 . 

In Table IVTI1 we present the combined fit results of the proton Y\{Q 2 ) data (elastic con- 
tribution excluded) in APT, the SGF model and conventional PT in PDG and denominator 
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forms. One can see there is noticeable sensitivity of the extracted a 1 ™ and //4 with respect 
to the minimal fitting scale Q 2 min variations, which may be (at least, partially) compensated 
by their RG logQ 2 evolution, similar to the BSR case. For completeness we included in 
Table [VTT1 APT fits for a^ v (Ql) and ^(Qq), taking into account their RG evolution. 

TABLE VII: Combined fit results of the proton T^(Q 2 ) data (elastic contribution excluded). APT fit 
results arj and fif^g (at the scale Qq = lGeV 2 ) are given without and with taking into account the 
RG Q 2 evolution of a (Q 2 ) and fif PT (Q 2 ). 



Met nod 


Qmin> GeV 


a 


/ 71 /T2 


/ 71 /f 4 


/ 71 /rfi 

fi s /M° 




0.59 


0.33(3) 


-0.050(4) 








NLO PDG 


0.35 


0.43(5) 


-0.087(9) 


0.024(5) 







0.29 


0.37(5) 


-0.060(15) 


-0.001(8) 


0.006(5) 




0.59 


0.35(3) 


-0.058(4) 








NLO Denom 


0.20 


0.38(3) 


-0.076(4) 


0.013(1) 







0.17 


0.38(4) 


-0.070(8) 


0.010(4) 


0.0004(5) 




0.47 


0.32(4) 


-0.056(4) 








NLO SGF 


0.17 


0.36(3) 


-0.071(4) 


0.0082(9) 





M = 1 GeV 


0.10 


0.40(4) 


-0.080(4) 


0.0134(9) 


-0.0007(6) 




0.47 


0.35(4) 


-0.054(4) 








NNLO APT 


0.17 


0.39(3) 


-0.069(4) 


0.0081(8) 





no evolution 


0.10 


0.43(3) 


-0.078(4) 


0.0132(9) 


-0.0007(5) 




0.47 


0.33(4) 


-0.051(4) 








NNLO APT 


0.17 


0.31(3) 


-0.059(4) 


0.0098(8) 





with evolution 


0.10 


0.32(4) 


-0.065(4) 


0.0146(9) 


-0.0006(5) 



As we already mentioned, the evolution of the fJ^(Q 2 ) is taken to be the same as for the 
nonsinglet term ^^(Q 2 ), allowing one to keep only one fitting parameter /^(Qo) instead 
of two in the general case. We also tested that the singlet anomalous dimension instead of 
the nonsinglet one [resulting in the same Q 2 evolution of tfl(Q 2 ) as that of fJ^ +n (Q 2 )} leads 
to close fit results within error bars. 

Figure [H] demonstrates the characteristic values of the proton data fits x 2 jD.o.f. (upper 
row) and the twist-4 coefficient /i4 (lower row) as functions of ao at different values of 
(numbers at the curves). One can see that at lower Q 2 {Q 2 min < 1 GeV 2 ) the APT description 
(left panels) turns out to be more precise and stable than that in the standard PT (right 
panels). Though we have taken the fitted values of ao and higher twists [i2% in the minima of 
each x 2 /D.o.f. curve as best fits, the naive constraint x 2 /D.o.f. < 1 (dotted horizontal lines 
mark 1) provides a quite wide spread in the allowable values of the fit parameters. However, 
it would be reasonable to take the spread between different minima as an optimistic error 
bar of our analysis. This gives us the following result: a = 0.33 ± 0.05, which is consistent 
with the recent analysis by COMPASS Q and HERMES H (see Table U). 

In Fig. [91 we show the best fit results for the less precise neutron r™(Q 2 ) data. Again, 
the APT fit gives the HT values demonstrating a better convergence of the OPE series, in 
contrast to the usual PT fit results. Fits with APT and more precise Denom PT couplings 
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FIG. 8: Behavior of x 2 /D.o.f. and from the proton data fits (with only one 1/Q 2 term) as 
functions of ao at different values of Q 2 min (the numbers at the curves) in the APT (left panels) 
and PT (right panels) cases. 



lead to a much smaller value of yJl and more stable fitting curves than that with the PDG 
coupling. Also the axial singlet charge ao extracted within APT from the neutron data turns 
out to be very close to the one extracted from more precise proton data (see Table rvTTj) . 
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Q, GeV 



FIG. 9: Best (2+l)-parametric fits of the JLab and SLAC data on T™ calculated with NLO Denom (solid 
line) and PDG (dashed line) couplings and NNLO APT (dash-dotted line). 

To obtain the genuine value of the twist-4 term /i^, we act in a similar way as for the 
BSR case in the previous subsection, namely, subtracting an extra (Q 2 )~ 1 term induced by 



16 



the APT series 



A (1) / i 



Ens m = E NS (a s = a^{Q<)) + < 



Q 2 + °{q^ 



Ej PT (Q 2 ) = E s (a s = + < 



A (1 f V 

LO(n2\\ i .JS eff 



Q 2 



(4.2) 



with Aty f ~ 0.18 GeV, xf s = 2.035, and xf = 0.661, and using the fit result in APT (with 



evolution) ^ /M 2 = —0.065 from Table IVTij we obtain 



/4(lGeV 2 
M 2 



w* {^ APT + U a3 + l as ) xf 5A ^ 2 + \ <v xf A ^ 2 1 " "°- 055 - (4 - 3) 



Analogously, for a neutron we have /i^/M 2 ~ —0.010. Subtracting it from the proton 
value (I4.3p . we get for the nonsinglet twist-4 term ^~ n /M 2 ~ —0.045 , which is close to the 
result in Ref. [38|, showing up the consistence of the APT analysis with the usual PT one for 
the proton and neutron SSR Y p { n data at Q 2 > 1 GeV 2 . Our result ( 14. 3 p is also consistent 
with the previous extraction at higher energies in Ref. [36| within the error bars (see also 
Table HI). 

It is worth noting that the best APT fit allows one to describe low-energy JLab data 
on r^ ,n at scales down to Q ~ 350 MeV with only the first three terms of the OPE series, 
unlike the usual PT case, where such fits happened to be impossible (due to the ghost issue) 
even for an increasing number of HT terms. This means that the lower bound of the pQCD 
applicability (supported by power HT terms) now may be shifted down to Q ~ Aqcd — 350 
MeV. 

However, it seems to be difficult to get a description in the region Q < Aq C o- This is 
not surprising, because the expansion in positive powers of Q 2 and its matching [9] with the 
HT expansion are relevant here. In this respect, the Aqcd scale appears as a natural border 
between "higher-twist" and "crural" nonperturbative physics. 
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PT 
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FIG. 10: Scale dependence of the gluon polarization Ag, obtained for different versions of perturbation 
theory - in APT (solid line), in conventional PT (dashed line), and in the SGF model (dash-dotted line). 



Finally, in Fig. [TUJ we show the scale dependence of the gluon polarization Ag obtained 
in APT, PT, and the SGF model. In conventional PT the value of Ag is small at the lower 
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scale Q 2 ~ 0.3 GeV 2 (see Ref. |55|). However, as one can see from Fig. [10], one may evolve Ag 
starting from higher scales Q 2 > 1 GeV 2 down to the deep infrared region and observe that 
the smallness of Ag is a consequence of the Landau singularities in a s . Applying different 
ghost-free models we see that Ag is much higher at Q 2 < 0.5 GeV 2 than one predicted in 
the standard PT. 



V. CONCLUSION AND OUTLOOK 

The singlet axial charge ao is the essential element of the nucleon spin structure which 
is related to the average total quark polarization in the nucleon. In this paper, we system- 
atically extracted this quantity from very accurate JLab data on the first moments of spin 
structure functions g\ ,n . 

These data were obtained at low Q 2 region 0.05 < Q 2 < 3 GeV 2 , and therefore, a special 
attention was paid to the QCD coupling in this domain. We demonstrated that the denom- 
inator form (12.101) of the QCD coupling a s is more suitable at the low Q 2 (see Figs. [1] and 
[2]). In particular, at the two-loop level it happens to be quite close to the exact numerical 
solution of the corresponding two-loop RG equation for Q > 0.5 GeV. 

The performed analysis includes even lower Q ~ ^qcd and involves the QCD coupling 
which is free of Landau singularities. For this purpose we used the APT (TiJ and the soft 
glueball-freezing model [l7[ for the infrared-finite QCD coupling a s . It was shown that the 
singularity-free APT and SGF QCD couplings are very close in the domain Q > 400 MeV. 

One can argue that large order perturbative and nonperturbative contributions are mixed 
up, and the duality between them is expected (see Ref. j56|). We tested a separation of 
perturbative and nonperturbative physics and performed a systematic comparison of the 
extracted values of the higher-twist terms in different versions of perturbation theory. A 
kind of duality between higher orders of PT and HT terms is observed so that higher order 
terms absorb part of the HT contributions moving the pQCD frontier between the PT and 
HT contribution to lower Q values in both nonsinglet and singlet channels (see Fig. [3]). As 
expected, the value of a changes substantially when coming from LO to NLO, whereas it 
is quite stable in higher-loop approximations. 

The perturbative contribution to the proton spin sum rule and to the Bjorken sum 
rule r^ _n in the APT approach and the SGF model is less than 5 % for Q > A. This explains 
the similarity of the extracted higher-twist parameters for these two modifications of QCD 
couplings. 

In the APT approach the convergence of both the higher orders and HT series is much 
better. In both the nonsinglet and singlet case, while the twist-4 term happened to be 
larger in magnitude in the APT than in the conventional PT, the subsequent terms are 
essentially smaller and quickly decreasing (as the APT absorbs some part of nonperturbative 
dynamics described by HT). This is the main reason for the shift of the pQCD frontier 
to lower Q values. A satisfactory description of the proton SSR and BSR data down to 
Q ~ Aqcb ~ 350 MeV was achieved by taking the higher-twist and (analytic) higher 
order perturbative contributions into account simultaneously (see Figs. [6] and [7]). The best 
accuracy for the extracted values of ao and higher-twist contributions «2i is achieved for the 
most precise proton SSR data while the analysis of the data on the neutron SSR shows the 
compatibility with the analysis of the BSR which is free from the singlet contribution. 

For the first time we considered the QCD evolution at low Q 2 of both the leadings twist 
a and the higher-twist « 4 terms using the (fractional) analytic perturbation theory [54j and 
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also the related evolution of the average gluon polarization Ag. Account of this evolution, 
which is most important at low Q 2 , improves the stability of the extracted parameters 
whose Q 2 dependence diminishes (see Table IVlT]) . As a result, we extract the value of the 
singlet axial charge a (l GeV 2 ) = 0.33 ± 0.05. This value is very close to the corresponding 
COMPASS 0.35 ± 0.06 H and HERMES 0.35 ± 0.06 @ results. 

The RG evolution of ao is related to the evolution of the average gluon polarization Ag 
[l|, 0]. The results of the evolution of Ag in the analytic perturbation theory and in the 
standard PT was compared (see Fig. [TO]) . The decrease of Ag at low Q 2 in APT is not so 
dramatic as in the standard PT case 551 ] . 

In a sense, it could be natural if the main reason for the significant shift of the pQCD 
frontier to lower Q 2 scales was the disappearance of unphysical singularities in perturbative 
series. Note that the data at very low Q ~ Aqcd are usually dropped from the analysis of 
ao and the higher-twist term in the standard PT analysis because of Landau singularities. 
At the same time, the compatibility of our results for a , extracted from the low energy 

29| demonstrates the universality of the nucleon spin 



JLab data with previous results [2 
structure at large and low Q 2 scales. It will be very interesting to explore the interplay 
between perturbative and nonperturbative physics against other low energy experimental 
data. 
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